#include "stdinc.h"
#include "newhash.h"
#include "extfunc.h" 
#include "extvab.h" 

static ubyte binLight[16384];
static ubyte probableMatrix[4][4] = {
// A  C  T  G								 7 3 2 1
	{7, 2, 1, 3},				//	A->A G C T
	{2, 7, 3, 1},				//	C->C T A G
	{1, 3, 7, 2},				//	T->T C G A
	{3, 1, 2, 7}				//	G->G A T C	
};

/* static ubyte2 doubleBitMasker[7] = { */
/* 	0x3,								//000000 00000011 */
/* 	0xC,								//000000 00001100 */
/* 	0x30,								//000000 00110000 */
/* 	0xC0,								//000000 11000000 */
/* 	0x300,							//000011 00000000 */
/* 	0xC00,							//001100 00000000 */
/* 	0x3000							//110000 00000000 */
/* }; */

/* static boolean staticFlag=1; */

static long long readsInGap=0;

static int buffer_size=10000000;
static long long readCounter;
static long long mapCounter;
long long single_count;
long long single_map;
static int ALIGNLEN=0;
//buffer related varibles for chop kmer
static int read_c;
static char **rcSeq;
static char **seqBuffer;
static int *lenBuffer;
static unsigned int *ctgIdArray;
static int *posArray;
static char *orienArray;
static char *footprint;   // flag indicates whether the read shoulld leave markers on contigs

// kmer related variables
static int kmer_c;
static Kmer *kmerBuffer,*hashBanBuffer;
static kmer_t **nodeBuffer;
static boolean *smallerBuffer;
static unsigned int *indexArray;

static int *deletion;

static void parse1read(int t,int threadID);
static void threadRoutine(void *thrdID);
static void searchKmer(int t,KmerSet *kset);
static void chopKmer4read(int t,int threadID);
static void thread_wait(pthread_t *threads);

static void creatThrds(pthread_t *threads,PARAMETER *paras)
{
	unsigned char i;
	int temp;
	
	for(i=0;i<thrd_num;i++){
		//printf("to create %dth thread\n",(*(char *)&(threadID[i])));
		if((temp = pthread_create(&threads[i], NULL, (void *)threadRoutine, &(paras[i]))) != 0){ 
			printf("create threads failed\n");
			exit(1);
		}
	}
	//printf("%d thread created\n",thrd_num);
}

static void threadRoutine(void *para)
{
	PARAMETER *prm;
	int i,t;
	unsigned char id;
	
	prm = (PARAMETER *)para;
	id = prm->threadID;
	//printf("%dth thread with task %d, hash_table %p\n",id,prm.task,prm.hash_table);
	while(1){
		if(*(prm->selfSignal)==1){
			for(i=0;i<kmer_c;i++)
			{
				//if((hashBanBuffer[i]&taskMask)!=prm.threadID)
				if((hashBanBuffer[i]%thrd_num)!=id)
					continue;
				searchKmer(i,KmerSets[id]);
			}
			*(prm->selfSignal) = 0;
		}else if(*(prm->selfSignal)==2){
			for(i=0;i<read_c;i++){
				if(i%thrd_num!=id)
					continue;
				chopKmer4read(i,id+1);
			}
			*(prm->selfSignal) = 0;
		}else if(*(prm->selfSignal)==3){
		// parse reads
			for(t=0;t<read_c;t++){
				if(t%thrd_num!=id)
					continue;
				parse1read(t,id+1);
			}
			*(prm->selfSignal) = 0;
		}else if(*(prm->selfSignal)==5){
			*(prm->selfSignal) = 0;
			break;
		}

                usleep(1);
	}
}
/*
static void chopReads()
{
	int i;
	for(i=0;i<read_c;i++){
		chopKmer4read(i,0);
	}
}
*/
static void chopKmer4read(int t,int threadID)
{
	int len_seq=lenBuffer[t];

	if(len_seq<overlaplen+1)
		return;
	
	char *src_seq=seqBuffer[t];
	char *bal_seq=rcSeq[threadID];
	int j,bal_j;
	Kmer hash_ban,bal_hash_ban;
	Kmer word,bal_word;
	int index;
	
	word = 0;
	for (index = 0;index<overlaplen;index++){
		word <<= 2;
		word += src_seq[index];
	}
	reverseComplementSeq(src_seq, len_seq,bal_seq);

		// complementary node
	bal_word = reverseComplement(word,overlaplen);
	bal_j = len_seq-0-overlaplen;  //  0;     
	index = indexArray[t];
	if(word<bal_word){
		hash_ban = hash_kmer(word);
		kmerBuffer[index] = word;
		smallerBuffer[index] = 1;
		hashBanBuffer[index++] = hash_ban;
	}else{
		bal_hash_ban = hash_kmer(bal_word);
		kmerBuffer[index] = bal_word;
		smallerBuffer[index] = 0;
		hashBanBuffer[index++] = bal_hash_ban;
	}

		//printf("%dth: %p with %p\n",kmer_c-1,bal_word,bal_hash_ban);
	for(j = 1; j <= len_seq - overlaplen; j ++)	{
		word = nextKmer(word,src_seq[j-1+overlaplen]);
		bal_j = len_seq-j-overlaplen; //  j; 
		bal_word = prevKmer(bal_word,bal_seq[bal_j]);
		if(word<bal_word){
			hash_ban = hash_kmer(word);
			kmerBuffer[index] = word;
			smallerBuffer[index] = 1;
			hashBanBuffer[index++] = hash_ban;
				//printf("%dth: %p with %p\n",kmer_c-1,word,hashBanBuffer[kmer_c-1]);
		}else{
				// complementary node
			bal_hash_ban = hash_kmer(bal_word); 
			kmerBuffer[index] = bal_word;
			smallerBuffer[index] = 0;
			hashBanBuffer[index++] = bal_hash_ban;
				//printf("%dth: %p with %p\n",kmer_c-1,bal_word,hashBanBuffer[kmer_c-1]);
		}
	}
}

//splay for one kmer in buffer and save the node to nodeBuffer
static void searchKmer(int t,KmerSet *kset)
{
	kmer_t *node;
	
	boolean found = search_kmerset(kset, kmerBuffer[t], &node);

	++kset->searchCnt;
	
	if(found)
	{
		++kset->foundCnt;
		if(!node->deleted)
			nodeBuffer[t] = node;
		else
		{
			++kset->delCnt;
			nodeBuffer[t] = NULL;
		}
	}
	else
	{
		++kset->searchSpcSeedCnt;

		boolean spcFlag;
		Kmer buff_kmer, spc_kmer;
		ubyte2 spc_bases;
		spcKmer *rs;
		spcBase *tmpBase;
		
		buff_kmer=kmerBuffer[t];
		spc_kmer = ((buff_kmer>>14)&0xFFFFFFF00) | ((buff_kmer>>12)&0xC0) | ((buff_kmer>>10)&0x3C) | ((buff_kmer>>6)&0x3);
		spc_bases = ((buff_kmer>>8)&0x3000) | ((buff_kmer>>6)&0xC00) | ((buff_kmer>>2)&0x3C0) | (buff_kmer&0x3F);
		
		spcFlag = search_spckmerset(spcSet, spc_kmer, &rs);
		
		if(spcFlag)
		{
			++kset->getSpcSeedCnt;
			
			int j=0,getFlag=-1;
			int mismatch=0;
			ubyte2 tmp,mostLastBase = 0;	//loci flags
			ubyte2 bestSpcBases = 0;			//best spaced bases
			int min_mis=31;
			ubyte2 tmpSpcBase;
							
			tmpBase=rs->start;
			
			//fprintf(stderr,"search %llu\tspc_kmer %u\tspc_bases %u\n", kmerBuffer[t], spc_kmer, spc_bases);
			
			while(tmpBase != NULL)
			{
				tmpSpcBase = tmpBase->spaced_bases;
				tmp = ((spc_bases ^ tmpSpcBase) & 0x5555) | (((spc_bases ^ tmpSpcBase) & 0xAAAA)>>1);
				mismatch=binLight[tmp];
				
				if(mismatch < min_mis)				//get the minimal mismatch spaced_bases
				{
					min_mis = mismatch;
					mostLastBase = tmp;
					bestSpcBases = tmpSpcBase;
					node = tmpBase->large_kmer;
					getFlag=0;
				}
				else if(mismatch == min_mis)	//if same amount of mismatch, choose the most right mismatch pos
				{
					if(tmp<mostLastBase)
					{
						mostLastBase = tmp;
						bestSpcBases = tmpSpcBase;
						node = tmpBase->large_kmer;
						getFlag=1;
					}
					else if(tmp == mostLastBase)	//if same mismatch pos, choose the most probable one[see probableMatrix]
					{
/*
static ubyte probableMatrix[4][4] = {
//A  C  T  G							 7 3 2 1
7, 2, 1, 3,				//	A->A G C T
2, 7, 3, 1,				//	C->C T A G
1, 3, 7, 2,				//	T->T C G A
3, 1, 2, 7				//	G->G A T C	
};
*/
						getFlag=2;
						ubyte2 readBases = spc_bases, loopBases = tmpSpcBase, bestBases = bestSpcBases, mismatchFlag = tmp;
						for(j=0;j<7;j++)
						{
							if((mismatchFlag & 0x3) > 0)
							{
								if(probableMatrix[(readBases & 0x3)][(loopBases & 0x3)] > probableMatrix[(readBases & 0x3)][(bestBases & 0x3)])
									//check each 2 bits(1 base) if mismatch
								{
									mostLastBase = tmp;
									bestSpcBases = tmpSpcBase;
									node = tmpBase->large_kmer;
									break;
								}
								else if((probableMatrix[(readBases & 0x3)][(loopBases & 0x3)] < probableMatrix[(readBases & 0x3)][(bestBases & 0x3)]))
									break;
							}
							mismatchFlag>>=2;
							readBases>>=2;
							loopBases>>=2;
							bestBases>>=2;
						}
					}
				}
				
				tmpBase = tmpBase->next;
			}

			if(getFlag<0)
			{
				fprintf(stderr,"getFlag error at %llu",kmerBuffer[t]);
				exit(-1);
			}
			++kset->levelGet[getFlag];
			nodeBuffer[t] = node;
		}
		else
			nodeBuffer[t] = NULL;
	}
}

static void parse1read(int t,int threadID)
{
	unsigned int j,i,s;
	unsigned int contigID;
	int counter2=0,counter;
	unsigned int ctgLen,pos;
	kmer_t *node;
	boolean isSmaller;
	int flag,maxOcc=0;
	int alldgnLen = lenBuffer[t] > ALIGNLEN ? ALIGNLEN:lenBuffer[t];
	int multi = alldgnLen-overlaplen+1 < 5 ? 5:alldgnLen-overlaplen+1;
	unsigned int start,finish;

	footprint[t] = 0;

	start = indexArray[t];
	finish = indexArray[t+1];
	if(finish==start){  //too short
		ctgIdArray[t] = 0;
		deletion[threadID]++;
		return;
	}
	for(j=start;j<finish;j++)
		if(nodeBuffer[j])
			counter2++;
	if(counter2<2)
		deletion[threadID]++;
	counter = counter2 = 0;
	for(j=start;j<finish;j++){
		node = nodeBuffer[j];	
		if(!node)  //same as previous
			continue;
		flag = 1;
		for(s=j+1;s<finish;s++){
			if(!nodeBuffer[s])
				continue;
			if(nodeBuffer[s]->l_links==node->l_links){
				flag++;
				nodeBuffer[s] = NULL;
			}
		}
		if(flag>=2)
			counter2++; //a loose alignment 
		if(flag>=multi)
			counter++;
		else
			continue;
		if(flag>maxOcc){
			pos = j;
			maxOcc = flag;
		}
	}
	if(!counter){  //no match
		ctgIdArray[t] = 0;
		return;
	}
	if(counter2>1)
		footprint[t] = 1;  //aligned to multi contigs

	j = pos;
	i = pos - start + 1;
	node = nodeBuffer[j];	
	isSmaller = smallerBuffer[j];
	contigID = node->l_links;
	ctgLen = contig_array[contigID].length;
	pos = node->r_links;
	if(node->twin==isSmaller){
		orienArray[t] = '-';
		ctgIdArray[t] = getTwinCtg(contigID);
		posArray[t] = ctgLen - pos -overlaplen -i + 1;
	}else{
		orienArray[t] = '+';
		ctgIdArray[t] = contigID;
		posArray[t] = pos - i + 1;
	}
	
}

static void sendWorkSignal(unsigned char SIG,unsigned char *thrdSignals)
{
	int t;

	for(t=0;t<thrd_num;t++)
		thrdSignals[t+1] = SIG;
	while(1){
		usleep(10);
				
		for(t=0;t<thrd_num;t++)
			if(thrdSignals[t+1])
				break;
		if(t==thrd_num)
			break;
	}
}

/* static void locate1read(int t) */
/* { */
/* 	int i,j,start,finish; */
/* 	kmer_t *node; */
/* 	unsigned int contigID; */
/* 	int pos,ctgLen; */
/* 	boolean isSmaller; */

/* 	start = indexArray[t]; */
/* 	finish = indexArray[t+1]; */
		
/* 	for(j=start;j<finish;j++){ */
/* 		node = nodeBuffer[j];	 */
/* 		if(!node)  //same as previous */
/* 			continue; */
/* 		i = j - start + 1; */
/* 		isSmaller = smallerBuffer[j]; */
/* 		contigID = node->l_links; */
/* 		ctgLen = contig_array[contigID].length; */
/* 		pos = node->r_links; */
/* 		if(node->twin==isSmaller){ */
/* 			ctgIdArray[t] = getTwinCtg(contigID); */
/* 			posArray[t] = ctgLen - pos -overlaplen -i + 1; */
/* 		}else{ */
/* 			ctgIdArray[t] = contigID; */
/* 			posArray[t] = pos - i + 1; */
/* 		} */
/* 	} */
	
/* } */

static void output1read(int t, FILE *outfp)
{
	int len = lenBuffer[t];
	int index;
	readsInGap++;
/*
	if(ctgIdArray[t]==735||ctgIdArray[t]==getTwinCtg(735)){
		printf("%d\t%d\t%d\t",t+1,ctgIdArray[t],posArray[t]);
		int j;
		for(j=0;j<len;j++)
			printf("%c",int2base((int)seqBuffer[t][j]));
		printf("\n");
	}
*/
	for(index=0;index<len;index++)
		writeChar2tightString(seqBuffer[t][index],rcSeq[1],index);
	fwrite(&len,sizeof(int),1,outfp);
	fwrite(&ctgIdArray[t],sizeof(int),1,outfp);
	fwrite(&posArray[t],sizeof(int),1,outfp);
	fwrite(rcSeq[1],sizeof(char),len/4+1,outfp);

}

static void getReadIngap(int t,int insSize,FILE *outfp,boolean readOne)
{
	int read1,read2;
	if(readOne){
		read1 = t;
		read2 = t + 1;
		ctgIdArray[read1] = ctgIdArray[read2];
		posArray[read1] = posArray[read2] + insSize - lenBuffer[read1];  //   --> R2       <-- R1
		output1read(read1,outfp);
	}else{
		read2 = t;
		read1 = t - 1;
		ctgIdArray[read2] = ctgIdArray[read1];
		posArray[read2] = posArray[read1] + insSize - lenBuffer[read2];  // --> R1     <-- R2
		output1read(read2,outfp);
	}
}

static void recordLongRead(FILE *outfp)
{
	int t;

	for(t=0;t<read_c;t++){
		readCounter++;
		if(footprint[t])
			output1read(t,outfp);
	}
}

static void recordAlldgn(FILE *outfp,int insSize,FILE *outfp2)
{
	int t,ctgId;
	boolean rd1gap,rd2gap;

	for(t=0;t<read_c;t++){
		readCounter++;
		single_count++;
		rd1gap = rd2gap = 0;
		ctgId = ctgIdArray[t];
		if(outfp2&&t%2==1){ //make sure this is read2 in a pair
			if(ctgIdArray[t]<1&&ctgIdArray[t-1]>0){
				getReadIngap(t,insSize,outfp2,0);  //read 2 in gap
				rd2gap = 1;
			}
			else if(ctgIdArray[t]>0&&ctgIdArray[t-1]<1){
				getReadIngap(t-1,insSize,outfp2,1);  //read 1 in gap
				rd1gap = 1;
			}
		}
		if(ctgId<1)
			continue;
		mapCounter++;
		single_map++;
		fprintf(outfp,"%lld\t%u\t%d\t%c\n",readCounter,
						ctgIdArray[t],posArray[t],orienArray[t]);
		if(t%2==0)
			continue;
		if(outfp2&&footprint[t-1]&&!rd1gap)
			output1read(t-1,outfp2);
		if(outfp2&&footprint[t]&&!rd2gap)
			output1read(t,outfp2);
		
	}
}

//load contig index and length
void basicContigInfo(char *infile)
{
	char name[256],lldne[1024];
	FILE *fp;
	int length,bal_ed,num_all,num_long,index;
	
	sprintf(name,"%s.ContigIndex",infile);
	fp = ckopen(name,"r");

	fgets(lldne,sizeof(lldne),fp);
	sscanf(lldne+8,"%d %d",&num_all,&num_long);
	//printf("%d edges in graph\n",num_all);
	num_ctg = num_all;
	contig_array = (CONTIG *)ckalloc((num_all+1)*sizeof(CONTIG));

	fgets(lldne,sizeof(lldne),fp);
	num_long = 0;
	while(fgets(lldne,sizeof(lldne),fp)!=NULL){
		sscanf(lldne,"%d %d %d",&index,&length,&bal_ed);
		
		contig_array[++num_long].length = length;
		contig_array[num_long].bal_edge = bal_ed+1;
		if(index!=num_long)
			printf("basicContigInfo: %d vs %d\n",index,num_long);
		if(bal_ed==0)
			continue;
		contig_array[++num_long].length = length;
		contig_array[num_long].bal_edge = -bal_ed+1;
			
	}

	fclose(fp);
}

void prlRead2Ctg(char *libfile,char *outfile)
{
	long long i;
	char *src_name,*next_name,name[256];
	FILE *fo,*outfp2=NULL;
	int maxReadNum,libNo,prevLibNo,insSize = 0;
	boolean flag,pairs=1;
	pthread_t threads[thrd_num];
	unsigned char thrdSignal[thrd_num+1];
	PARAMETER paras[thrd_num];

	maxReadLen = 0;
	maxNameLen = 256;
	scan_libInfo(libfile);
	alloc_pe_mem(num_libs);
	if(!maxReadLen)
		maxReadLen = 100;
	//printf("In file: %s, max seq len %d, max name len %d\n\n",
		//libfile,maxReadLen,maxNameLen);
	if(maxReadLen>maxReadLen4all)
		maxReadLen4all = maxReadLen;

//////////////////////////////////////////// spcSet
	fflush(stdout);

	ubyte2 spc_i,spc_j;
	for(spc_i=0;spc_i<16384;spc_i++)
	{
		binLight[spc_i]=0;
		for(spc_j=spc_i;spc_j;spc_j=spc_j&(spc_j-1))
			++binLight[spc_i];
	}

	spcSet = init_spckmerset(KmerSets[thrd_num-1]->size*thrd_num, 0.77f);

	
	for(i=0;i<thrd_num;i++)
	{
		buildSpcKmerSet(KmerSets[i], spcSet);
		mvnv(0,"%lldth spaced bases set build complete.",i);
	}
//////////////////////////////////////////// END spcSet

	src_name = (char *)ckalloc((maxNameLen+1)*sizeof(char));
	next_name = (char *)ckalloc((maxNameLen+1)*sizeof(char));

	kmerBuffer = (Kmer *)ckalloc(buffer_size*sizeof(Kmer));
	hashBanBuffer = (Kmer *)ckalloc(buffer_size*sizeof(Kmer));
	nodeBuffer = (kmer_t **)ckalloc(buffer_size*sizeof(kmer_t *));
	smallerBuffer = (boolean *)ckalloc(buffer_size*sizeof(boolean));

	maxReadNum = buffer_size/(maxReadLen-overlaplen+1);
	maxReadNum = maxReadNum%2==0 ? maxReadNum:maxReadNum-1;   //make sure paired reads are processed at the same batch
	seqBuffer = (char **)ckalloc(maxReadNum*sizeof(char*));
	lenBuffer = (int *)ckalloc(maxReadNum*sizeof(int));
	indexArray = (unsigned int *)ckalloc((maxReadNum+1)*sizeof(unsigned int));
	ctgIdArray = (unsigned int *)ckalloc((maxReadNum+1)*sizeof(unsigned int));
	posArray = (int *)ckalloc((maxReadNum+1)*sizeof(int));
	orienArray = (char *)ckalloc((maxReadNum+1)*sizeof(char));
	footprint = (char *)ckalloc((maxReadNum+1)*sizeof(char));

	for(i=0;i<maxReadNum;i++)
		seqBuffer[i] = (char *)ckalloc(maxReadLen*sizeof(char));

	rcSeq = (char **)ckalloc((thrd_num+1)*sizeof(char*));
	deletion = (int *)ckalloc((thrd_num+1)*sizeof(int));
	thrdSignal[0] = 0;
	deletion[0] = 0;
	if(1){
		for(i=0;i<thrd_num;i++){
			rcSeq[i+1] = (char *)ckalloc(maxReadLen*sizeof(char)); 
			deletion[i+1] = 0;
			thrdSignal[i+1] = 0;
			paras[i].threadID = i;
			paras[i].mainSignal = &thrdSignal[0];
			paras[i].selfSignal = &thrdSignal[i+1];
		}
		creatThrds(threads,paras);
	}
	
	if(!contig_array)
		basicContigInfo(outfile);
	
	sprintf(name,"%s.readInGap",outfile);
	outfp2 = ckopen(name, "wb");

	sprintf(name,"%s.readOnContig",outfile);
	fo = ckopen(name,"w");
	fprintf(fo,"read\tcontig\tpos\n");
	readCounter = mapCounter = 0;
	single_map = single_count = 0;
	kmer_c = n_solexa = read_c = i = libNo = readNumBack = gradsCounter = readsInGap = 0;
	prevLibNo = -1;
	
	//mvnv(0,"Start loading reads.");
	
	while((flag=read1seqInLib(seqBuffer[read_c],next_name,&(lenBuffer[read_c]),&libNo,pairs,0))!=0)	{
		if(libNo!=prevLibNo){
			prevLibNo = libNo;
			insSize = lib_array[libNo].avg_ins;
			ALIGNLEN = lib_array[libNo].map_len;
			
			if(insSize>1000)
				ALIGNLEN = ALIGNLEN < 35 ? 35: ALIGNLEN;
			else
				ALIGNLEN = ALIGNLEN < 32 ? 32: ALIGNLEN;
			//printf("current insert size %d, map_len %d\n",insSize,ALIGNLEN);
			
		}
	
		if(insSize>1000)  
			ALIGNLEN = ALIGNLEN < (lenBuffer[read_c]/2+1) ? (lenBuffer[read_c]/2+1):ALIGNLEN;
		
//		if((++i)%100000000==0)
//			printf("[%s]%lld reads processed.\n",__FUNCTION__,i);
		indexArray[read_c] = kmer_c;
		if(lenBuffer[read_c] >= overlaplen+1)     
			kmer_c += lenBuffer[read_c] - overlaplen + 1;
		read_c++;
		if(read_c==maxReadNum){
			//mvnv(0,"Start processing reads.");
			
			indexArray[read_c] = kmer_c;
			
			sendWorkSignal(2,thrdSignal);
			//mvnv(0,"chop finished one buffer.");
			sendWorkSignal(1,thrdSignal);
			//mvnv(0,"search finished one buffer.");
			sendWorkSignal(3,thrdSignal);
			//mvnv(0,"parse finished one buffer.");

			recordAlldgn(fo,insSize,outfp2);
			kmer_c = 0;
			read_c = 0;
		}
	}

	if(read_c){
		indexArray[read_c] = kmer_c;
		sendWorkSignal(2,thrdSignal);
		sendWorkSignal(1,thrdSignal);
		sendWorkSignal(3,thrdSignal);
		recordAlldgn(fo,insSize,outfp2);
		//printf("Output %lld out of %lld (%.1f)%% reads in gaps\n",readsInGap,readCounter,
			//		(float)readsInGap/readCounter*100);
	}
	if(readCounter)
		printf("[%s]total %llu reads , map-rate (%.1f)%%\n",__FUNCTION__,
			readCounter,(float)mapCounter/readCounter*100);
	sendWorkSignal(5,thrdSignal);
	
	thread_wait(threads);
	fclose(fo);

	sprintf(name,"%s.peGrads",outfile);	
	fo = ckopen(name,"w");
	fprintf(fo,"grads&num: %d\t%lld\t%d\n",gradsCounter,n_solexa,maxReadLen4all);
	if(pairs){ 
		if(gradsCounter)
			;
			//printf("%d pe insert size, the largest boundary is %lld\n\n",
				//gradsCounter,pes[gradsCounter-1].PE_bound);
		else
			printf("no paired reads found\n");
		for(i=0;i<gradsCounter;i++)
			fprintf(fo,"%d\t%lld\t%d\t%d\n",pes[i].insertS,pes[i].PE_bound,pes[i].rank,pes[i].pair_num_cut);
	}
	fclose(fo);

	fclose(outfp2);
	
	free_pe_mem();
	free_libs();
	if(1){ // multi-threads
		for(i=0;i<thrd_num;i++){
			deletion[0] += deletion[i+1];
			free((void *)rcSeq[i+1]);
		}
	}
	//printf("%d reads deleted\n",deletion[0]);
	
	ubyte8 searchCntTot=0, foundCntTot=0, delCntTot=0,
					searchSpcSeedCntTot=0,getSpcSeedCntTot=0,levelGet1=0,levelGet2=0,levelGet3=0;
	for(i=0;i<thrd_num;i++){
		searchCntTot += KmerSets[i]->searchCnt;
		foundCntTot += KmerSets[i]->foundCnt;
		delCntTot += KmerSets[i]->delCnt;
		searchSpcSeedCntTot += KmerSets[i]->searchSpcSeedCnt;
		getSpcSeedCntTot += KmerSets[i]->getSpcSeedCnt;
		levelGet1 += KmerSets[i]->levelGet[0];
		levelGet2 += KmerSets[i]->levelGet[1];
		levelGet3 += KmerSets[i]->levelGet[2];		
	}
	fprintf(stderr,"SEARCH: Search %llu, get %llu, deleted %llu\n", 
									searchCntTot, foundCntTot, delCntTot);	
	fprintf(stderr,"SPACED SEED: Search %llu, get %llu, LVnum %llu, LVpos %llu, LVpro %llu\n", 
									searchSpcSeedCntTot, getSpcSeedCntTot, levelGet1, levelGet2, levelGet3);
	
	free((void *)rcSeq);
	free((void *)deletion);
	for(i=0;i<maxReadNum;i++)
		free((void *)seqBuffer[i]);
	free((void *)seqBuffer);
	free((void *)lenBuffer);
	free((void *)indexArray);
	
	free((void *)kmerBuffer);
	free((void *)smallerBuffer);
	free((void *)hashBanBuffer);
	free((void *)nodeBuffer); 
	free((void *)ctgIdArray);
	free((void *)posArray);
	free((void *)orienArray);
	free((void *)footprint);

	free((void *)src_name);
	free((void *)next_name);

	if(contig_array){
		free((void *)contig_array);
		contig_array = NULL;
	}
}

static void thread_wait(pthread_t *threads)
{
	int i;
	for(i=0;i<thrd_num;i++)
		if(threads[i]!=0)
			pthread_join(threads[i],NULL);
	
}

/********************* map long reads for gap filling ************************/
void prlLongRead2Ctg(char *libfile,char *outfile)
{
	long long i;
	char *src_name,*next_name,name[256];
	FILE *outfp2;
	int maxReadNum,libNo,prevLibNo;
	boolean flag,pairs=0;
	pthread_t threads[thrd_num];
	unsigned char thrdSignal[thrd_num+1];
	PARAMETER paras[thrd_num];

	maxReadLen = 0;
	maxNameLen = 256;
	scan_libInfo(libfile);
	if(!maxReadLen)
		maxReadLen = 100;
	int longReadLen = getMaxLongReadLen(num_libs);
	if(longReadLen<1)  // no long reads
		return;
	maxReadLen4all = maxReadLen<longReadLen ? longReadLen:maxReadLen;
	//printf("In file: %s, long read len %d, max name len %d\n\n",
		//libfile,longReadLen,maxNameLen);
	maxReadLen = longReadLen;

	src_name = (char *)ckalloc((maxNameLen+1)*sizeof(char));
	next_name = (char *)ckalloc((maxNameLen+1)*sizeof(char));

	kmerBuffer = (Kmer *)ckalloc(buffer_size*sizeof(Kmer));
	hashBanBuffer = (Kmer *)ckalloc(buffer_size*sizeof(Kmer));
	nodeBuffer = (kmer_t **)ckalloc(buffer_size*sizeof(kmer_t *));
	smallerBuffer = (boolean *)ckalloc(buffer_size*sizeof(boolean));

	maxReadNum = buffer_size/(maxReadLen-overlaplen+1);
	maxReadNum = maxReadNum%2==0 ? maxReadNum:maxReadNum-1;   //make sure paired reads are processed at the same batch
	seqBuffer = (char **)ckalloc(maxReadNum*sizeof(char*));
	lenBuffer = (int *)ckalloc(maxReadNum*sizeof(int));
	indexArray = (unsigned int *)ckalloc((maxReadNum+1)*sizeof(unsigned int));
	ctgIdArray = (unsigned int *)ckalloc((maxReadNum+1)*sizeof(unsigned int));
	posArray = (int *)ckalloc((maxReadNum+1)*sizeof(int));
	orienArray = (char *)ckalloc((maxReadNum+1)*sizeof(char));
	footprint = (char *)ckalloc((maxReadNum+1)*sizeof(char));

	for(i=0;i<maxReadNum;i++)
		seqBuffer[i] = (char *)ckalloc(maxReadLen*sizeof(char));

	rcSeq = (char **)ckalloc((thrd_num+1)*sizeof(char*));
	deletion = (int *)ckalloc((thrd_num+1)*sizeof(int));
	thrdSignal[0] = 0;
	deletion[0] = 0;
	if(1){
		for(i=0;i<thrd_num;i++){
			rcSeq[i+1] = (char *)ckalloc(maxReadLen*sizeof(char)); 
			deletion[i+1] = 0;
			thrdSignal[i+1] = 0;
			paras[i].threadID = i;
			paras[i].mainSignal = &thrdSignal[0];
			paras[i].selfSignal = &thrdSignal[i+1];
		}
		creatThrds(threads,paras);
	}
	
	if(!contig_array)
		basicContigInfo(outfile);
	
	sprintf(name,"%s.longReadInGap",outfile);
	outfp2 = ckopen(name, "wb");

	readCounter = 0;
	kmer_c = n_solexa = read_c = i = libNo;
	prevLibNo = -1;
	while((flag=read1seqInLib(seqBuffer[read_c],next_name,&(lenBuffer[read_c]),&libNo,pairs,4))!=0)	{
		if(libNo!=prevLibNo){
			prevLibNo = libNo;
			ALIGNLEN = lib_array[libNo].map_len;
			ALIGNLEN = ALIGNLEN < 35 ? 35: ALIGNLEN;
			//printf("Map_len %d\n",ALIGNLEN);
		}
	
//		if((++i)%100000000==0)
	//		printf("%lld reads processed.\n",i);
		indexArray[read_c] = kmer_c;
		if(lenBuffer[read_c] >= overlaplen+1)     
			kmer_c += lenBuffer[read_c] - overlaplen + 1;
		read_c++;
		if(read_c==maxReadNum){
			indexArray[read_c] = kmer_c;
			
			sendWorkSignal(2,thrdSignal);
			sendWorkSignal(1,thrdSignal);
			sendWorkSignal(3,thrdSignal);

			recordLongRead(outfp2);
			kmer_c = 0;
			read_c = 0;
		}
	}

	if(read_c){
		indexArray[read_c] = kmer_c;
		sendWorkSignal(2,thrdSignal);
		sendWorkSignal(1,thrdSignal);
		sendWorkSignal(3,thrdSignal);
		recordLongRead(outfp2);
		//printf("Output %lld out of %lld (%.1f)%% reads in gaps\n",readsInGap,readCounter,
			//		(float)readsInGap/readCounter*100);
	}

	sendWorkSignal(5,thrdSignal);
	
	thread_wait(threads);

	fclose(outfp2);
	
	free_libs();
	if(1){ // multi-threads
		for(i=0;i<thrd_num;i++){
			deletion[0] += deletion[i+1];
			free((void *)rcSeq[i+1]);
		}
	}
	//printf("%d reads deleted\n",deletion[0]);
	
	free((void *)rcSeq);
	free((void *)deletion);
	for(i=0;i<maxReadNum;i++)
		free((void *)seqBuffer[i]);
	free((void *)seqBuffer);
	free((void *)lenBuffer);
	free((void *)indexArray);
	
	free((void *)kmerBuffer);
	free((void *)smallerBuffer);
	free((void *)hashBanBuffer);
	free((void *)nodeBuffer); 
	free((void *)ctgIdArray);
	free((void *)posArray);
	free((void *)orienArray);
	free((void *)footprint);

	free((void *)src_name);
	free((void *)next_name);

}

